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Abstract 

The spectral properties of the Wilson-Dirac operator in 2-dimensional QED 
responsible for the appearance of exceptional configurations in quenched sim- 
ulations are studied in detail. The mass singularity structure of the quenched 
functional integral is shown to be extremely compicated, with multiple branch 
points and cuts. The connection of lattice topological charge and exactly real 
eigenmodes is explored using cooling techniques. The lattice volume and spac- 
ing dependence of these modes is studied, as is the effect of clover improvement 
of the action. A recently proposed modified quenched approximation is applied 
to the study of meson correlators, and the results compared with both naive 
quenched and full dynamical calculations of the same quantity. 
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1 Introduction 



A basic quantity in essentially all lattice calculations of hadronic properties is the 
quark propagator in the presence of a fixed external gauge field. Hadron propagation 
amplitudes or current matrix elements are extracted from averages of products of 
quark propagators over ensembles of gauge fields generated relative to either a pure 
gluonic weight e'-^"^'''' ^^^^'^ action) ^j^g quenched approximation, or e'-^"^'''^ ^^^^'^ action) . 
(determinant of Dirac operator) for a fully dynamical calculation including internal 
quark loop effects. There has naturally been considerable attention devoted to the 
issue of the extent to which the quenched approximation deviates from the full the- 
ory [|l|. Are the differences of a detailed quantitative character only, perhaps mainly 
absorbable in a redefinition of the scale of the theory, or are there serious qualita- 
tive differences which make it difficult to extract reliable physics from a quenched 
calculation? 

Qualitative differences which distinguish the quenched theory from the full dy- 
namics concern the sensitive dependence of certain amplitudes on the quark mass m 
in the chiral limit m — 0. Generic effects include the "hairpin" diagrams which signal 
the generation of the singlet rj' mass and the corresponding "quenched chiral logs" 
which are expected to give logarithmic corrections to the mass dependence in the 
chiral limit. In lattice calculations, finite volume effects and convergence problems in 
evaluating the light quark propagators are also elements which impact the study of 
the chiral limit. A particular artifact associated with quenched calculations involving 
light Wilson fermions is the occasional appearance of "exceptional" gauge configura- 
tions with apparently nonsensical values for the light quark propagators resulting in 
extremely noisy statistical averages and unreliable values for the hadronic amplitudes. 
Clearly, it will be impossible to determine, in any definitive way, the chiral structure 
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of the quenched theory until the nature of the exceptional configuration problem is 
fully understood and the resulting pathologies tamed. Our objective here is to study, 
in a specific model - 2 dimensional QED- under very precise numerical control, the 
role of the exceptional configuration phenomenon in quenched calculations with light 
Wilson fermions. In the final section, we shall use this model to perform a series of 
numerical tests on a recently proposed cure of the exceptional configuration problem 
P which suggests that greatly improved results are indeed possible for calculations 
using Wilson-Dirac fermions in the low mass chiral regime. 

Following the example of Smit and Vink we can take advantage of the fact that 
2-dimensional QED (massive Schwinger model) provides a particularly useful testbed 
for investigating the pathologies of the quenched approximation. As in 4-dimensional 
QCD, there is an intimate connection between the topological structure of the theory 
and the spectrum of the Dirac operator near the origin, which controls the chiral limit 
behavior. Moreover, proposed modifications of the quenched approximation can be 
tested against the full theory {including determinant effects) in great detail in the 
2-dimensional case. 

On large lattices at weak coupling the Wilson spectrum {A^} settles (taking QED2 
as an example) near the real axis into four distinct branches at Re(A) = ±2, 0, with 
the central region containing two branches. It is usual to select the undoubled left 
branch (see Fig(l)) as a model for the continuum theory by tuning the bare quark 
parameter m appropriately. For example, in QED2 one finds a critical line at 
Re(A) ^ — 2 + 0.65//9, with a dispersion of the real part of the eigenvalue around 
this value which grows rapidly as one goes to stronger coupling. In particular, at 
stronger coupling one finds an increasing density of exactly real eigenmodes straying 
to the left away from the critical line into the region corresponding to small positive 
quark masses. Consequently, evaluation of the propagator in configurations giving 
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rise to such modes necessarily results in noisy results, as one occasionally ends up 
sitting very close to a pole of the propagator. In a sense, the original rationale for 
moving from Minkowski to Euclidean space - avoiding the poles of the propagator 
at Minkowski momenta- has been vitiated by the complex (i.e. non-skewhermitian) 
character of the Wilson-Dirac operator. The exceptional configurations observed in 
quenched calculations appear to be due to precisely this feature of the Wilson action, 
and the study of these configurations reduces essentially to a study of the frequency 
and location of the exactly real eigenvalues of the Wilson-Dirac matrix defining the 
quadratic fermion action of a lattice theory. 

The present paper is organized as follows. In Section 2 we review some of the 
basic properties of the Wilson-Dirac matrix in lattice gauge theory. Most of the re- 
sults listed here are well known but the review allows us to clarify our notation and 
to remind the reader of some exact properties of the Wilson-Dirac spectrum follow- 
ing purely from linear-algebraic considerations. Results are given here both for QED 
in 2 and QCD in 4 dimensions. Section 3 contains a more detailed analysis of the 
mass singularity structure of the quenched functional integral. It is shown that the 
quenched functional integral as defined in a Monte Carlo simulation is strictly speak- 
ing undefined. The functional integral must be defined by analytic continuation from 
a pole- free region of the complex quark mass plane. In the region of physical interest, 
the path integral defines a cut function in this plane with a complicated branch struc- 
ture. Section 4 describes the relationship between topological charge and the exactly 
real modes responsible for the pathologies of the quenched approximation. In Section 
5 we study the dependence of the real part of the Wilson-Dirac spectrum on lattice 
volume and lattice spacing. In Section 6 we give the results of a detailed study of the 
changes in the real spectrum induced by the addition of a clover improvement term 
IP to the action. Finally, in Section 7 we test the validity of a recently proposed 
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modified quenched approximation (MQA) which resolves the singularities introduced 
by real eigenmodes by a simple pole-shifting procedure. In particular, we carry out 
an explicit comparison of the MQA with both the naive quenched and full dynamical 
results for correlators of quark bilinears. 

2 Formal properties of the Wilson-Dirac operator 

The Euclidean Dirac operator I^{A) in the continuum is formally skew-hermit ian 
for smooth gauge fields A, so at a formal level the quark propagator {I^{A) + m)^^ 
exists for all nonzero quark mass m. The existence of the inverse at zero quark mass is 
related to the topological structure of the background gauge field by the Atiyah-Singer 
index theorem. For example, in 2D QED, defining the Dirac operator on a compact 
2-sphere of radius R by the usual one-point compactification |0| gives a rigorously 
self-adjoint iI^{A) (for gauge fields which approach pure gauge at Euclidean infinity) 
with a purely discrete spectrum : 

lI^{A)Xa{x) = XaXaix) (1) 
/ ^^^ ^2^j^2 ^i(^)^fe(^) = ^''b (2) 

e 2R 

-^e^^F^^ + ^2^ji2 xhbXaix) = pure divergence (3) 

The eigenfunctions Xi ci^re normalizable relative to the inner product (2), and the 
"pure divergence" referred to in (3) is just the divergence of the axial vector current 
of the model, which is seen to reduce to a sum of a topological term and a sum over 
the 75 eigenvalues of each of the discrete modes. The latter sum (as a consequence of 
{75, I^{A)} = 0) reduces to a sum over zero modes only, and therefore to a difference 
of the number of positive and negative chirality zero mode solutions of the Dirac 
equation in the background field. Nielsen and Schroer have given an explicit analytic 
expression for these solutions in an arbitrary smooth (and asymptotically pure-gauge) 



field. Essentially the same statements hold for QCD in 4 dimensions, with the obvious 
replacement of e^^F^j^y by the appropriate 4-dimensional Chern-Simons FF term. 

Once the Euclidean space-time is discretized as in lattice theory, little of the pre- 
ceding formalism remains in exact form, though qualitative remnants of the zero mode 
structure can be observed provided one works on large lattices with very smooth back- 
ground gauge fields. In the first place, although a discrete spectrum is now guaranteed 
(the Dirac operator is a finite matrix!), the addition of a Wilson term to eliminate 
doubling destroys the skew-hermitian nature of the Dirac operator. Specifically, quark 
propagators are inverses of a matrix D — rW + m = Ai + m, with D, W and m the 
naive Dirac matrix, W the Wilson term, and m a quark mass parameter: 



M = D-rW (4) 

Daarh,bl3n = fji) abU ap{rn^i)5n,rh+[L - ^i) abU^^^ (5) 

Waarh,bl3n = ^Sab{Ual3{'mfJ')Sfi,m+fi + Ulp{w)^H,m^ (6) 

Here a,b are Dirac indices, a, (3 color indices, m,n lattice sites, and r the Wilson 
parameter, which must be chosen < r < 1 and is usually taken to be unity. As D 
is skew-hermitian and W hermitian, the full Wilson-Dirac Ai = D — rW matrix is 
complex, and the generic eigenvalue (which locates poles of the quark propagator in 
the complex mass plane) will also be complex. 

It will frequently be convenient to use a single index i = {aarh) to identify the 
single Grassmann variable for a quark field component, and to write the Wilson- 
Dirac matrix elements with the notation A4ij. From (^ it follows immediately that 
although Ai is neither skewhermitian (unless r=0) nor hermitian, the matrix j^Ai is 
hermitian. Consequently the coefficients of the secular equation for Ai are real: 

det(A-75A<) = (det(A*-75-M))* (7) 
6 



det{X-M) - {det{X*-M)y (8) 

so that if A is an eigenvalue of A^, so is A*- the eigenvalues are therefore either real 
or appear in conjugate pairs. The nearest neighbour structure of Al also implies 
immediately that 

TriM"") = 0, n odd (9) 

which implies that odd powers of A are absent in the secular equation det(A — A4) = 0. 
Consequently, if A is an eigenvalue, so is —A. The generic case therefore has the roots 
of the secular equation for M. appearing as quartets a, —a, a*, —a* or as pairs of 
real eigenvalues r, — r. To summarize, we may write (the dimension is 12V for 4D 
QCD, 2V for 2D QED, where V is the lattice volume = number of lattice sites) 



det(A -M)=X^ + C2A^-2 + C4A^-^ + ...cn (10) 
where the Cj are all real and given in terms of Ai by standard trace formulas 

C2 = -lTr{M') (11) 

C4 = -l{Tr{M')-lTr{Mr) (12) 

C6 = -l{Tr{M^) + lTr{M^f-^Tr{M^)Tr{M')), etc (13) 
6 8 4 



Henceforth we suppose that r=l. The factors 1 ± 7^ appearing in M. are then 
projection operators, and in particular (1 + ^f^){l — 7^^) vanishes. Consequently, all 
terms in Tr(AI^), which necessarily correspond to products A4ijA4ji where i,j involve 
nearest neighbour lattice sites, automatically vanish. As a result 

C2 = Tr{M^) = (14) 

and the sum of squares of the complex eigenvalues of the Wilson-Dirac matrix must 
vanish in an arbitrary gauge configuration! This sum rule is very useful in checking 
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the accuracy of the (sometimes unstable) numerical routines we have used to effect 
the spectral resolution of A4. The calculation of C4 is equally straightforward. Again, 
each contributing term in Tr(A^'^) = A4ijA4jkAiki-Mii corrresponds to a sequence of 
4 links forming a closed path and with (for r=l) backtracks forbidden. The resulting 
color trace simply produces the plaquette trace for an elementary square on the lattice, 
provided that the linear dimension of the lattice (in all directions) is greater than 4, 
so that Polyakov lines stretching across the whole lattice are excluded. In D=2,4 
dimensions, one finds 

Tr(A^^) = -L»^ReTr([/p) (15) 

p 

C4 = ^Y.^eTr{Up) (16) 
^ p 

which gives the sum rule that the fourth power of eigenvalues of M. must sum to the 
gauge action. Similarly, C6,C8, .. involve sums of non-backtracking gauge-invariant 
loops (and eventually, Polyakov lines) of length 6,8, etc. Note the symmetry property 
under simultaneous inversion of all links (which will be used below in discussing the 
analytic structure of quenched amplitudes) Cn{{Ui}) — Cn{{Uj^^}). 

Similar calculations show that the quantity Tr(A^ A^^) is independent of the gauge 
field configuration. Namely 

Tt{MM^) = W (2D QED) (17) 
Tr(MM^) = 16A^c^ (4D QCD) (18) 

This implies, for example, that the spectrum of M. is rigorously contained within a 
circle of radius 2L for 2 dimensional lattice QED on a LxL lattice. A stronger bound 
follows from the observation that (for r=l) TrAI" < C - 3"^ for large n as the number 
of closed non-backtracking paths of length n is clearly < 3" in 2 dimensions and the 
Dirac traces involve products of projection operators which do not grow with n. This 
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Figure 1: A typical Wilson-Dirac spectrum in QED2 

implies |Aa| < 3 for eigenvalues ol Ai. In fact with typical gauge field configurations, 
eigenvalues are not found with magnitude larger than 2. A typical spectrum for 
quenched 2D QED (12x12 lattice at /3=4) is shown in Fig(l). The left, right and 
central branches of the Wilson-Dirac spectrum are clearly visible, as well as a pair of 
exactly real eigenvalues for each branch. One of these real modes has drifted quite 
a substantial distance into the interior of the spectral oval, away from the critical 
regime. 

The discrete symmetries of the Wilson matrix mean that the spectrum is consid- 
erably more constrained than would be the case for a generic complex operator. In 
particular the appearance of exactly real eigenvalues as in Fig(l) is not a miracle. 
As a simple toy example, consider A^=4, and assume that the secular equation is 
constrained, as for our Wilson matrix, to have only even terms with real coefficients. 
For a general value of the r parameter (not equal to unity), we would have C2 7^ 0. 
The resulting 4 roots are all real provided only that C2 < 0, > — 4c4 > 0. The 
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reality of such roots is clearly robust with respect to small variations of 02,04. (For 
r=l, similar arguments show the robustness of the reality of roots of the cubic sec- 
ular equation for A^=6). Explicit calculations (described below) of the full spectrum 
in 2D lattice QED on a variety of lattices confirm that when real eigenvalues of Ai 
appear for some set of link angles, they remain real in a finite neighborhood of this 
point in the configuration space of gauge link variables. The appearance of exactly 
real modes may seem unsurprising in view of the necessity of such modes in a topo- 
logically nontrivial background in the continuum, but such modes do not appear in 
the Kogut-Susskind (KS) formulation of staggered fermions on the lattice [§, where 
the fermion matrix (essentially D of (5)) is exactly skew-hermit ian, so that the spec- 
trum is confined to the imaginary axis, and an exactly zero eigenvalue is nongeneric. 
Instead, the latticized KS matrix displays pure imaginary eigenvalues close to, but 
not exactly, zero in a smooth background lattice field (e.g. after cooling) of nonzero 
winding number. 

3 Mass Singularity Structure of the Quenched Func- 
tional Integral 

The appearance of poles in the physical region for bare quark masses is a symptom 
of a serious disease in the standard quenched functional integral. Strictly speaking, 
this integral - even on a finite lattice, with compact integrations over gauge link 
variables- does not exist for values of the bare quark mass m in the physical range. 
The quenched functional integral can in principle be defined by analytic continuation 
from a pole-free region (see below) in the complex plane of m, but as we shall soon 
show the analytic structure of quenched amplitudes becomes extremely complicated 
in the region of physical interest. In any event analytic continuation is out of the 
question for a numerically computed amplitude. We are forced to face the unpleasant 
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fact that the quenched functional integral with Wilson fermions, as defined in a 
manner amenable to Monte Carlo sampling, runs directly over poles of the integrand. 
Staying at large quark masses would avoid this problem, which only appears, as we 
shall see, at kappa values k > Kq, where Kq is the critical kappa value for the free 
theory (=0.25 in 2 dimensions, 0.125 in 4 dimensions). On the other hand, Monte 
Carlo simulation of the quenched integral at smaller, physically interesting bare quark 
mass will necessarily give ill-defined results, even when very large ensembles are used. 

In a full dynamical simulation of hadron correlators in lattice QCD, one computes 
the functional integral over both fermionic and gauge degrees of freedom of a product 
of an even number of quark fields 



Such 2n-point functions are necessarily meromorphic functions of the bare quark 
mass provided the determinant arising from the Grassmann integration is properly 
included. Indeed, by the usual rules of Grassmann integration, on a finite lattice 
both the numerator N{m) and the vacuum amplitude Zf)['m) in (13) arc polynomials 
in the quark mass. In the quenched approximation, on the other hand, the absence 
of the determinantal factor results in a much more complicated analytic structure 
in the complex m-plane, as poles in quark propagators Dij = {A4 + m)^j^ remain 




uncancelled by the corresponding zeroes in det(A4 + m). The position of these poles 
in the complex m-plane is of course dependent on the gauge field and consequently 
the location of singularities in m of N^m), where the gauge field integration has 
been performed, depends on whether the gauge field integrals are prevented from 
deformation by either a pinch or an end point singularity. (Of course, in the quenched 
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approximation, the normalizing denominator ZQ{m) only involves the pure gauge 
action and is consequently independent of m, so the entire singularity structure derives 
from the numerator in (13).) The generic structure of a quenched contribution to 
N{m) is thus 

/ ^i^^ir^, (rT\ N-2 I \rr^ jV-4 , 777\^ ' (analytic in m, Ui) (21) 

where q is the number of quark propagators appearing in the correlator. Here we 
are examining the analytic structure in the complex plane of a single quark mass- 
the issue of several flavors, and the effect of hairpin graphs in isosinglet sectors in 
reducing the level of singularity in (|2T|) will be discussed below. As the coefficients 
Cm{Ui) are all gauge-invariant, a maximal tree of links may be assumed fixed to unity 
and the integrals in (|2T|) only include physical degrees of freedom. 

For the rest of this section we shall consider abelian gauge theory on a 2 dimen- 
sional LxL lattice, although much of the discussion applies just as well to 4D QCD. 
The gauge integrations are then over link angles < 6i < 27t,1 < p where p is the 
number of physical (non-fixed) links. From the bound on the magnitude of eigen- 
values of Ai implied by (p!7D, it is clear that N{m) is analytic for |m| > 2L, and 
that singularities are confined to a compact central region around the origin of the 
m plane with poles (if any) and paired branch points arranged symmetrically around 
the origin and real axis . Alternatively, defining zi = e*^' , the link integrations can be 
thought of as over a p-dimensional torus T, with singularities in m arising whenever 
two or more root surfaces + C2{zi)m^~'^ + ...cn{zi) = pinch T from opposite 
sides. That branch point singularities can easily arise in this fashion is readily seen 
from a simple example (involving only a single physical link) where the integral over 
the link variable z = e*^ is written as a complex integral over the unit circle 

/M^/t^ ^.. , , . = . .1 „ . (22) 



27Tizm^ + a{z + l/z)m'^ + 1 ^(l + m^)^ -Aa 
12 



Here a is a real constant. The branch point singularities located at the roots of 
(1 + m^)^ = Aa^m^ arise from a pinch of the unit circle integration contour by the 
two roots in z of the integrand in (0). More generally, the roots of det(m — M.{zi)) 
define singularity surfaces of the integral 



1 -r-r dZl 



— Uirh (24) 



This integral is well-defined for |m| > 2L, and defines a real-analytic function of m 
there, but will encounter difficulties as m is brought into the central region where roots 
of det(m — Ai{zi)) can occur. In particular, branch points of f{m) in the complex 
plane of the mass variable m will occur at points tuq whenever the hypertorus of 
integration in ( PB]) is pinched from opposite sides by two branches of the singularity 
surface Si{zi,mo) and 5*2 (-2/, mo) with 

ttiT^ + a2T^ = (25) 
dzi dzi 

for fixed ai,a2 and all links zi. For example, setting all zi = 1 (the ordered free 
configuration), there is (for periodic boundary conditions) always a pair of real roots 
at tuq = ±2. If mo is moved infinitesimally away from these points (say to mo = 2±e) 
roots of det(m — Ai{zi)) appear with \zi\ < 1 and \zi\ > 1, where the configuration 
{zf^} is necessarily a root if {zi} is (as Cn{zi) = Cn{zf'^)). Thus the singularity 
surfaces are naturally paired by taking S2{zi,m) = Si{zf^,m), in which case ( pSD 
is automatically satisfied with ai = 0:2 • The surfaces come together and pinch the 
integration hypertorus at a single point provided zi = zf^ for all links, i.e. if zi = ±1. 
We conclude that branch points of /(m) are present at all complex values of m which 
solve det(m — Ai{zi)) = for any configuration where the link variables are either 
1 or -1. In particular, the eigenvalues of Ai for the ordered configuration are all 
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Figure 2: Location of some of the branch points for L=6 

branch points of the quenched amphtudes. In addition, there are many other branch 
points corresponding to configurations where the hnk variables are allowed to take the 
values ±1. A glance at Fig(3), which shows branch points arising from the ordered 
configuration as well as from a disordered configuration where the link variables were 
set to ±1 randomly (on a 6x6 lattice), indicates just how complicated the analytic 
structure of a lattice quenched amplitude can be. In particular, branch cuts extend 
along the entire real axis from m = —2 to m = 2, exactly in the region where one 
attempts to approach continuum behavior. 
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4 Zero modes, topological charge, and cooling 

As discussed above in Section 2, there is a clear connection in the continuum theory 
between the topological charge of the background gauge configuration and the zero 
mode structure of the Dirac operator. On the lattice, this connection is necessarily 
less precise, as discretization effects smear out both the shape of the Dirac spectrum 
and alter the value of the topological charge associated with a given lattice gauge 
field. As discussed by Smit and Vink, two possible definitions of topological charge 
for 2D QED on the lattice are 

Qi = 7^Esin(M (26) 

where Op is the plaquette angle for plaquette P. This definition of course gives in 
general a noninteger value for the topological charge. Alternatively, one may define 

Qs^^^^P, -7r<^p<7r (27) 

which gives an integer value for the topological charge, agreeing with p6|) for smooth 
fields in the continuum limit. For a typical gauge configuration at weaker values of the 
coupling the number of exact zero modes found for the lattice Wilson-Dirac operator 
agrees with Q2 (in each critical branch), with Q2 essentially the closest integer to Qi, 
while at lower (3 one finds increasingly frequent discrepancies in these three measures 
of topological charge. 

A convenient way of understanding the role of discretization on a given finite lattice 
involves cooling a Monte Carlo generated configuration to gradually remove high- 
momentum components from the gauge configuration. Such cooling tends to stabilize 
any nontrivial topology while bringing the zero mode and topological properties of the 
system closer to continuum-like behavior. A typical result of a series of cooling sweeps 
(in which the action is adiabatically lowered, effectively by running the simulation 
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Table 1: Cooling a configuration with topological charge = 1 



Cooling sweeps 


Qi 


Q2 


Action 


Real modes (left branch) 


Chirality Fraction 





0.850 


1 


20.8 


-1.892 


0.930 


1 


0.913 


1 


14.2 


-1.901 


0.999 


2 


0.940 


1 


10.1 


-1.911 


0.999 


3 


0.950 


1 


6.9 


-1.919 


0.999 


4 


0.959 


1 


4.7 


-1.926 


0.999 



at a high fixed /9 value) on the near-real-axis spectrum of M. (left branch) is shown 
in Table(l). As the initial configuration, with topological charge (5i=0.85 is cooled, 
the action decreases steadily, the topological charge Q\ increases towards unity, and 
the norm fraction of the positive chirality part of the corresponding real eigenmode 
increases rapidly towards 1. There is a single exactly real eigenmode (in each branch) 
of the Wilson-Dirac matrix throughout. Restricting our attention to the left critical 
branch only, cooling moves the real mode gradually further to the left (as we would 
expect by forcing the system to a larger effective /5), but leaves the number of real 
modes unchanged. Note that the discrete definition Q2 gives exactly unity for the 
topological charge throughout. 

In the course of a quenched Monte Carlo simulation topological charge fluctuations 
arise resulting in the appearance of new exactly real eigenmodes of the Dirac operator. 
An interesting example of this phenomenon is illustrated in Fig. 3. Initially (circular 
points) there is one exactly real eigenvalue per branch, and a pair of slightly off-real- 
axis complex modes considerably to the right of the left-branch critical line, which 
for the lattice in question (12x12 at /3=5) is at Re(m)Ri —1.91. Both Qi and Q2 give 
values for this initial Monte Carlo configuration close to -|-1 {Q2 is of course exactly 
1), as shown in Table(2). The corresponding eigenvector would have exactly positive 
chirality in the continuum- instead the norm fraction of the positive chirality piece is 
initially about 98.8% (Table(2), column 6). As the cooling proceeds, the complex pair 
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Table 2: Effect of cooling on a topological charge fluctuation 



Cooling sweeps 


Qi 


Q2 


Action 


Real modes (left branch) 


Chirality Fraction 





1.167 


1 


42.4 


-1.90 


0.988 


2 


1.266 


1 


25.9 


-1.92 


0.994 


4 


1.486 


2 


15.7 


-1.94,-1.71 


0.995,0.462 


6 


1.636 


2 


10.3 


-1.95,-1.82 


0.995,0.673 


8 


1.693 


2 


7.45 


-1.96,-1.86 


0.995,0.771 


10 


1.817 


2 


5.33 


-1.96,-1.90 


0.997,0.897 



in Fig(3) move onto the real axis (diamond points in Fig(3)) and then split laterally, 
with one member moving out towards the critical line. The topological charge Q\ 
moves smoothly up towards 2 as this happens, with the discrete version Qi switching 
from 1 to 2 at roughly the point when the second exactly real mode appears (we 
take snapshots of the entire spectrum after each cooling sweep through the whole 
lattice, so it is difficult to be more precise as to the exact point at which the new 
real modes appear) . Note that the "right chirahty" fraction of the new zero mode is 
initially rather low (less than one-half), but increases steadily as this mode becomes 
stabilized and moves leftward into the critical region. 

Evidently, in this case the initial configuration, with Qi, Qi of order unity, actu- 
ally contains a "nascent" topological charge 2 structure, with large high momentum 
fluctuations hiding the larger value of topological charge. The presence of this hidden 
topological charge is nevertheless signalled by the presence of complex pairs in the 
interior of the spectral oval, close to the critical region. In the case illustrated here, 
the effect of cooling mimics the further evolution under Monte Carlo updating, which 
also leads to Qi ~ 2 after a few more sweeps, with a spectrum containing two exactly 
real modes per branch. In other cases, we have observed a reduction of topological 
charge under cooling: a pair of real eigenvalues in the left and central branches move 
out into the interior of the spectral oval, colliding there, then moving vertically away 
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Effect of cooling on zero mode structure 
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Figure 3: Effect of cooling on near-zeromodes (/3=5, 12x12 lattice) 



from the real axis, eventually to join the noncritical modes far from the real axis. 
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5 Lattice Volume and Lattice Spacing Dependence 
of Real Spectrum 



It has become part of the lore of quenched lattice QCD simulations that the ex- 
ceptional configuration problem appears to be ameliorated at fixed f3 (i.e. lattice 
spacing) when one goes to larger physical volumes. It is therefore of some interest to 
see whether this feature can be understood in terms of the volume dependence of the 
real part of the Wilson-Dirac spectrum. As one increases the physical volume, one 
naturally expects more frequent occurrence of configurations with nonvanishing topo- 
logical charge, so it is clear that the overall number of exactly real zero modes, which 
as we have seen above are clearly correlated with nontrivial topological structures, 
should increase with volume. The appearance of exceptional configurations requires 
however that these modes appear in the quark mass regime being studied, just to the 
left of the critical line corresponding to zero renormalized mass. 

The simulations described in the preceding section show that the number of exact 
zero modes is roughly correlated with the absolute magnitude of the topological charge 
Qi = ^Spsin(^p). Quenched compact 2 dimensional QED is an exactly solvable 
model in the pure gauge sector, so it is a straightforward exercise to evaluate the 
distribution of this quantity over quenched configurations generated on a LxL lattice 
(periodic boundary conditions) at any given /3. If the lattice volume is sufficiently 
large that {Ii{(3) / Iq{(3))^^ (which decreases exponentially fast with lattice volume) is 
negligible, where /q, h are the usual modified Bessel functions, then periodic and free 
boundary conditions become equivalent and the distribution of Q\ becomes exactly 
Gaussian. (For the smallest lattice used below, L=6 at /5=4, this approximation 
incurs an error of order 10~^). One finds 




(28) 
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which imphes that the average value < \Qi\ > rises hnearly with the lattice size, or 
as the square root of the lattice volume. We have evaluated the spectra for 2500 
configurations of lattices at P =4 and sizes L=6,8, and 10. The total number of 
exactly real modes rises in this range of lattice sizes somewhat faster than L though 
not as rapidly as - specifically, the average number of exactly zero modes per branch 
was found to be 0.26, 0.47 and 0.63 for the L=6,8, and 10 lattices respectively. The 
spectral histograms for these three cases in the region of the left critical branch is 
shown in Fig(4). 

250 T 1 




Re(m) 



Figure 4: Volume dependence of spectrum (left branch)) 

Fig(4) provides a possible clue to the observation that exceptional configurations 
become rarer at fixed quark mass as the volume increases. We see that although the 
total number of real modes is certainly increasing with lattice volume, the leftmost 
zero mode shifts towards the boundary of the spectral region at m = —2 rather slowly 
(for L—6,8, and 10, the leftmost modes were found at -1.920,-1.923, and -1.927 respec- 
tively) , while the peak of the histogram, which we may identify with the critical mass 
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Figure 5: Beta dependence of spectral histogram 

value, shifts leftward more rapidly (for L=6,8 and 10, the histograms peak at -1.855, 
-1.87, -1.88), so that the probability of encountering an exceptional configuration as 
one approaches the critical line from the left decreases with increasing volume if one 
keeps a fixed offset from the critical line to maintain a fixed physical quark mass. 
However, these histograms also make it clear that the exceptional configurations will 
necessarily appear at any volume once one goes sufficiently close to the critical point. 
As we shall now see, the benefits to be obtained from increasing lattice volume are far 
less dramatic than from the reduction of lattice spacing (at fixed physical volume). 

In Fig(5) we show a histogram of exactly real eigenmodes appearing in quenched 
simulations in QED2 on LxL lattices with L=6,8,10 and 12, this time with the gauge 
coupling readjusted to keep the physical lattice volume fixed (specifically, we have 
used P—1.6, 2.9, 4.5 and 6.5). The histograms were obtained from 1000 configu- 
rations separated by 200 quenched Monte Carlo sweeps (except for L=12, where 
measurements were separated by 500 sweeps) . Note that the direction of the horizon- 
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tal axis in these plots is reversed from the previous spectral plots, as by convention 
positive quark mass corresponds to the region to the left of the critical line (at roughly 
2-0.65//5). The density of real modes near the left critical branch is shown, relative to 
a critical line defined in this case by the peak of the distribution (this turns out to be 
close to but slightly offset from the Smit-Vink critical line value, with the offset de- 
creasing with increasing (3), and with the scale on the horizontal axis rescaled to keep 
physical masses fixed. It is apparent that the spread of dangerous poles into the phys- 
ical mass region becomes acute at strong coupling, and that (at least in 2D QED) the 
probability of encountering exceptional configurations at a fixed renormalized quark 
mass decreases rapidly as /5 is increased. 

The manifest asymmetry of the real mode distribution around the critical line 
accounts for the uncontrollably noisy results obtained if one attempts a quenched 
simulation at lighter quark masses approaching from the interior of the spectral oval 
(i.e. at values of the kappa parameter greater than kappa critical). Nevertheless, we 
shall see below that the modified quenched approximation approach to resolving 
the exceptional configurations in the k, < Kc regime is equally successful in removing 
this noise and restoring sensible correlators for k > k^. 
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Dispersion of real poles for clover-improved actions 




Real pole position 



Figure 6: Dependence of spectral histogram on clover improvement coefficient {(3 =4.5) 

6 Effects of clover improvement on the Wilson- 
Dirac Spectrum 

The dispersion of real axis eigenvalues of the Wilson-Dirac matrix was seen above to 
decrease rapidly with increasing /3, and it is clear that we are dealing with a pathology 
which is directly traceable to the particular discretization of the theory implied by an 
unimproved Wilson action for the quark fields. At first sight this would suggest that 
the problem might be removable (or at least, substantially reduced) by adopting a 
locally improved action, say by adding a clover-leaf term with appropriate coefficient. 
Unfortunately (see also [Q|) this does not seem to be the case. In Fig. 6 we show 
the effect of introducing such a term, with various values of the coefficient, on the 
spectral histogram of real eigenvalues for a 8x8 lattice at /3=4.5 in QED2. Evidently, 
the clover improvement term is simply not effective in restraining the tendency of real 
modes to drift away from the critical region at strong coupling. 
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Of course, the addition of a clover term certainly will shift the location of real 
modes on each individual configuration. Such a shift can have a dramatic effect on 
a given exceptional configuration, as it may be enough to move the propagator pole 
far enough from the chosen quark mass value to remove the large fluctuation in the 
correlators for that configuration. But in a statistical sense the problem remains. Our 
QED2 studies have revealed a rather interesting effect at the individual configuration 
level: it appears that the effect of a clover term is much stronger on the two central 
branches {not used in typical simulations) of the Wilson-Dirac spectrum. This effect 
is shown in Fig(7), where the full spectrum on a typical topological charge Q2 —2 
configuration is shown for clover coefficient values of 0,1.2, and 2.4 (where the tree 
value is unity). At higher values of the clover coefficient the real modes in the right 
central branch become degenerate and the nonzero modes line up quite accurately 
with the ones on the real axis. There is also a tendency for the two central branches to 
separate as the clover coefficient is made larger. The effect on the left and right lateral 
branches is less pronounced- the two zero modes on the left branch move somewhat 
apart, while those on the right move together. The whole process is reversed left to 
right if the sign of the clover coefficient is changed. 
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Figure 7: Spectrum of topological charge 2 configuration, clover coefficients =0,1.2,2.4 
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7 Testing the modified quenched approximation 

Recently, we have proposed a nonlocal improvement procedure which should become 
arbitrarily precise in the continuum limit, as it involves explicitly shifting the position 
of real poles of the quark propagator by amounts which go to zero in this limit, but 
in such a way as to preserve the location of the critical point to leading order as 
the lattice spacing goes to zero. A redefinition of the Wilson-Dirac matrix which 
involves a shift in the location of real eigenvalues while leaving the corresponding 
eigenmodes unchanged evidently corresponds to a nonlocal addition to the action, 
as the typical eigenmodes for a real eigenvalue will be associated with a topological 
charge fluctuation which extends over several lattice sites (with power rather than 
exponential fallofF away from the "instanton center"). Of course, we need to be sure 
that this redefinition does not introduce an unwanted sizeable systematic distortion 
of the quenched averages at finite lattice spacing- in particular, does the proposed 
redefinition, which is really a prerequisite to having any precise notion of a quenched 
lattice amplitude, actually yield results which are less noisy and closer to the results 
of a full dynamical simulation on the same lattice? Direct evidence for this will now 
be presented. 

The modified quenched approximation (MQA) recently proposed by us yields a 
well-defined regularization of the quenched functional integral which is constructed 
to reduce to the usual version in the continuum limit, in the sense that real axis poles 
of the Wilson-Dirac propagator which have strayed into the physical region for the 
quark masses are shifted back to the critical line where they would sit in this limit. 
The shift is performed in such a way as to preserve the location of the critical line 
in an average sense (for more details, see [§]). In the context of QED2, we define 
the MQA modified propagator as follows. The fermion propagator S = + 
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(with the bare quark mass defined as mo = — in terms of the usual kappa 
parameters, and with ~2.0-0.65//3 = Ac in QED2) has the spectral decomposition 

= V ^"f , (29) 

where Aq are the eigenvalues of A4, and the indices i,j incorporate spatial and Dirac 
information. The left-shifted real eigenvalues will cause poles at values of mo corre- 
sponding to K, < Kc- For QED2 we are able to perform a full spectral resolution for 
each configuration . Define u = 2mo — ~ ^ and replace any pole appearing at a 
position u — Up^ie ^ follows 

(30) 



U - Upole U U + Upoie 

At large mass (or large u) the first two terms in the expansion in 1/u are identical 
and terms linear in the shifts should average to zero. This procedure preserves the 
location of Kc to first order in the dispersion of real eigenmodes around the critical line 
(which is of course going to zero in the continuum limit). The full MQA propagator 
may be simply computed by adding a term to the naive fermion propagator Sij which 
incorporates the pole shift. Namely, for any configuration where a pole(s) appears in 
the physical region (operationally we have defined this in this paper to be the region 
between and 2mo on the real bare quark mass axis), we replace 

S^"^"- = Sij + AP°^^(«)ResP/'^ (31) 

with 

^P°le(^) = 1 1 1 (32) 

U U + Upole U - Upole 

Although the procedure outlined here is certainly valid in the continuum limit, 
where we expect the frequency of exceptional configurations at any fixed physical 
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quark mass to vanish, it is perhaps not apparent that the modification induced by 
the pole shift procedure does not introduce some large and misleading distortion of 
the correlators at stronger coupling. The level of numerical control possible with 
QED2 allows a direct comparison of the naive quenched, shifted pole, and full dy- 
namical results for the interesting correlators of the theory. We have measured the 
pseudoscalar correlator < ■?/'(x)75'?/'(a;)?/'(0)75'?/'(0) > at fairly strong coupling (3 =4.5 
on a 10x10 lattice for these three situations for a variety of quark masses close to 
the critical point of the left critical branch. The full dynamical simulations were 
done using an exact update algorithm |]TD[ where the determinant and fermion prop- 
agator are known exactly at every stage. Typically, measurements were performed 
after every sweep for 400 sweeps, and the statistical error computed by measuring an 
autocorrelation time for each correlator (the naive quenched correlators decorrelate 
very rapidly so autocorrelation times are on the order of 1 sweep in the unshifted 
quenched case, while for the shifted pole and dynamical correlators, the autocorre- 
lation times range from 3-10 sweeps). QED2 is a superrenormalizable theory with a 
finite bare coupling which fixes the mass scale of the theory, so a direct comparison 
of quenched and dynamical results can be made in this theory without the need for 
a compensating scale change. As we work on fairly small lattices, it is difficult to 
find mass plateaus suitable for extracting a meson mass, and we have decided to 
compare directly the pseudoscalar correlators, which are after all physical quantities 
with unambiguous values (for given (3 and lattice size) in the full unquenched theory. 

The comparison of naive quenched, MQA and full dynamical simulations for quark 
mass values 0.08, 0.10 and 0.115 are shown in Fig(8). The quark mass values quoted 
here indicate the offset from a critical value on a 10x10 lattice at (3=4.5 defined as the 
average real part of the left branch real eigenmodes with —2 < Re(A) < —1.5 . The 
increase in the naive quenched errors at smaller quark mass can be directly traced 
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to the increasing frequency of nearby real poles. In all cases, the MQA correlator 
is less noisy and closer to the full dynamical result. Statistical errors in the latter 
are related to autocorrelation problems, not the intrinsic noisiness of the data. Of 
course, for very small quark masses, we begin to see a deviation between MQA and 
full dynamical results, as the MQA does not contain the suppression of nontrivial 
topologies implicit in the zero mass full dynamical theory. 

Amazingly, the MQA procedure makes it possible (to our knowledge, for the first 
time) to obtain sensible results by approaching the critical line from the interior of the 
spectral oval, corresponding to negative quark mass with the conventional definition 
(i.e. the "supercritical" regime with k > k^). In the continuum, a negative quark 
mass may be removed by a 75 redefinition of the quark field, but on the lattice this 
region contains many more real modes extending farther away from the critical line 
(see Figs. 4,5,6), leading to extremely noisy quenched correlators (although of course 
the full dynamical theory is perfectly well-defined here). The plots shown in Fig(9) 
for the two point function of ■075V' show that the MQA is effective in this case also, 
although the increased density of real modes in this regime makes the MQA less 
effective at extremely small quark masses- the results shown are for quark masses of 
-0.15, -0.16 and -0.185 (as before, on a 10x10 lattice at /3=4.5), with the negative 
sign indicating that we are now working in the interior of the spectral oval . For 
small quark masses, the results in the naive quenched model are again dominated by 
noise, while the MQA and full dynamical results are close , with the statistical noise 
under control in the MQA even for quite small negative quark masses. For larger 
quark masses, where sensible results can be extracted from an uncorrected quenched 
calculation, the MQA correlators lie between the naive quenched and full dynamical 
results, as was the case previously in the subcritical domain. 
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Figure 8: Pseudoscalar correlators in QED2- quenched, MQA and full dynamical 
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8 Conclusions 



The analytic arguments and numerical results presented here provide strong evidence 
for the interpretation of exceptional configurations as an artifact of quenched Wilson 
lattice theory arising from the structure of the spectrum of the discrete Wilson-Dirac 
matrix, and in particular from the appearance of exact real eigenmodes for partic- 
ular gauge configurations close to the critical line corresponding to the chiral limit. 
Two dimensional QED turns out to be a particularly useful testbed for studying this 
problem for two reasons. First, the zero mode and topological charge structure of the 
theory closely mimics that of 4 dimensional QCD, at least insofar as those aspects ger- 
mane to the real eigenmode problem are concerned. Secondly, extensive simulations 
can be performed at comparatively trivial computational cost (all of the numerical 
calculations reported in this paper were performed on a Pentium PC), yielding com- 
plete information about all aspects of the spectral structure. It is also possible to 
perform full dynamical simulations (with exact evaluation of the determinant) in this 
theory to test the validity of proposed modifications of the quenched approximation. 

Much of the material presented in Sections 2 and 3 extends in a direct and obvious 
fashion to 4 dimensional QCD. There also, discrete symmetries of the Wilson-Dirac 
secular equation lead to the appearance of exactly real eigenmodes. The quenched 
functional integral as conventionally defined will be a cut function of the bare quark 
mass with an exceedingly complex branch structure. Pinch arguments along the lines 
of Section 3 can again be used to identify a subset of branch points. Unfortunately, 
the very useful and detailed information on the lattice spacing, lattice volume and 
clover coefficient dependence of the real part of the Wilson-Dirac spectrum (sum- 
marized in the histograms of Figs. 4,5, and 6) is not yet available for 4 dimensional 
QCD. However the results presented here for QED2 confirm in all essential details the 
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picture of the zero mode problem in QCD presented in In particular, the direct 
comparisons of meson correlators in the conventional quenched, modified quenched 
(MQA), and exact dynamical simulations are a gratifying confirmation of the effi- 
cacy of the MQA method in resolving the exceptional configuration problem at finite 
lattice spacing. 
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